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ABSTRACT 


Determining target bearing based on a passive acoustic signal typically relies on 
beamforming the signals from an array of sound pressure sensors. A major 
drawback, however, is the proportional increase in array aperture when dealing 
with low frequencies, such as the lengthy towed arrays used for anti-submarine 
warfare. This thesis demonstrates the use of a single acoustic vector sensor 
(Microflown Ultimate Sound Probe (USP)) to derive the target bearing by 
processing both the pressure as well as particle velocity information of an 
acoustic wave. 

Field experiments were set up to track commercial aircraft during their 
final approach before landing. Despite healthy signal-to-noise (SNR) ratios, 
significant challenges were faced in accurate real-time tracking. Post-processing 
frequently achieved better results, but required the beamformer to process a 
broader range of frequencies (typically 300-1000 Hz), instead of focusing on 
narrowband energy peaks. This was attributed to the effects of noise and bottom 
reflections (mainly from the concrete ground), as implied by the distinct Lloyd’s 
mirror patterns in the spectrograms. Notwithstanding, additional information such 
as target altitude and horizontal distance at the closest point of approach (CPA) 
could be determined from analyzing these patterns. 
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I. INTRODUCTION 


The use of acoustics in target detection is of interest to many fields and 
has conventionally been achieved via phase beamforming spatially distributed 
sensors. The sensors first capture the pressure signals of the incoming acoustic 
waves. The signals are then combined in a beamformer to determine the 
direction of arrival (DOA). Such methods, however, have a major drawback. In 
order to achieve good bearing resolution, the aperture of an array must be large 
compared with the wavelength. Furthermore, a single linear array has complete 
azimuthal symmetry making it impossible to distinguish which side of the array 
the target is on. In low-frequency applications, they usually result in large 
physical footprints with consequently high cost, limited mobility, and counter¬ 
detection issues. 

Acoustic vector sensors, on the other hand, have the potential of 
achieving better DOA estimates with smaller arrays. Each acoustic wave gives 
rise to pressure as well as particle velocity changes in the transmission fluid. 
Knowledge of both would allow the DOA to be determined, as shown by Hawkes 
(1998). The Microflown Ultimate Sound Probe (USP) sensor combines three 
orthogonal particle velocity sensors and one pressure microphone. This enables 
wideband detection of both the polar and azimuthal angle of the acoustic DOA 
with just a single sensor with certain limitations which will be discussed. 

Caulk (2009) demonstrated this using a hybrid array of two Microflown 
sensors and one pressure microphone. He successfully determined the DOA of 
two sound sources operating at two different frequencies in the anechoic 
chamber to within about 5° accuracy. However, experiments on contacts of 
interest under field conditions failed due to wind noise issues. 

This paper builds upon Caulk’s thesis to demonstrate the Microflown 
USP’s ability to track actual planes during takeoff and landing at an airport. 
Advances include the direct extraction of sensor data into Matlab for processing 

1 



(without having to use LabView), and a near real-time display of DOA as visual 
feedback. A foam wind shield is also utilized to overcome the wind noise that 
prevented the USPs from tracking actual targets previously. 

The rest of this paper is organized as follows: 

• Chapter II: Theoretical background to vector sensors and their 
beamforming. 

• Chapter III: Hardware and software used, details of apparatus and 
experimental setup. 

• Chapter IV: Signal processing algorithms and considerations. 

• Chapter V: Results and analysis. 

• Chapter VI: Conclusion and recommendations. 
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II. THEORY 


A. AXES SETUP 

Figure 1 shows a close-up of the Microflown USP acoustic vector sensor, 
together with the axes configuration for this thesis. The three particle velocity 
sensors (BLUE, GREEN, RED) are orthogonally placed to detect sound from the 
corresponding Cartesian axes, (j) and 0 are the azimuthal and polar angles 
respectively. 



Figure 1. Schematic close-up of Microflown USP and axes setup 

B. BEAMFORMING VECTOR SENSORS 

Each vector sensor includes four sensors, one for pressure and three for 
particle velocity: 

e, is the pressure sensor output 

e 2 is the BLUE velocity sensor output (z-axis) 
e 3 is the GREEN velocity sensor output (x-axis) 
e 4 is the RED velocity sensor output (y-axis) 
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Although there are always slight deviations from orthogonality amongst 
the three velocity channels, previous work showed that these deviations were 
small compared to the uncertainty in the beamformer output. Let (0 2 ,(j> 2 ) be the 
direction of the Maximum Response Axis (MRA) for the BLUE vector sensor, 
(# 3 ,^ 3 ) be that for GREEN, and (9 4 ,(f> 4 ) be that for RED. The unit vectors 

corresponding to the MRA of each of the sensors would then be given in terms of 
the direction cosines as: 

Equation Section 2 

u e2 = sin 0 2 cos <j) 2 x + sin 0 2 sin <j) 2 y + cos 0 2 z = z ( 2 . 1 ) 

u e2 = sin #3 cos+ sinsin^ v + cos^z = x (2.2) 

u e4 = sin 0 4 cos <j) 4 x + sin 0 4 sin <j) 4 y + cos 0 4 z = y (2.3) 

It should be noted that various beamforming methods exist, and this thesis 
retains the method used by Caulk. Prior to beamforming, the outputs of the 
velocity sensors must be adjusted to be commensurate with the amplitude and 
phase of the pressure sensor output. This is done by measuring the transfer 
functions between the pressure and particle velocity outputs. While this can 
theoretically be derived via the supplied calibration data of the Microflown USPs, 
it may fail to capture effects of our mounting apparatus, or other sources of 
interference or noise from the rest of the electronics. A comprehensive method 
is to conduct in-situ calibration with the sensor in its experimental setup, which 
will be discussed in the following chapter. 

Once the outputs of the individual sensors have been properly calibrated, 
they can be put through the beamforming algorithm. For a steer direction of 
( 0 s ,fi s ), the unit vector in the steer direction is given in terms of the direction 
cosines as: 

u s = sin 0 s cos (j) s x + sin Q s sin <f> s y + cos 0 s z (2.4) 

Each particle velocity output would then be weighted accordingly using 

this steering vector by taking their dot products: 
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n = 2,3,4 


(2.5) 


w„ = ll, ll„ 


The pressure output is omnidirectional, hence: 


u s m el = 1 


( 2 . 6 ) 


Next, the velocity channels need to be expressed in terms of their 
equivalents to the pressure signal. This can be achieved by applying the 
appropriate transfer functions, described more in the next section. Eventually 
when given Xc n (k) : the discrete Fourier transform of ^corrected to the 
equivalent of the pressure sensor, the linear (Bartlett) beamformer output is given 
by: 

= ( 27 ) 


B(k,6 s ,</> s ) = 




1 = 1 


( 2 . 8 ) 


For a chosen k th frequency bin, the corresponding steer vector u s that 
generates the largest B(k,O s ,</> s ) would give the DOA of the incoming signal. 

Note (2.8) is only valid for a single vector sensor. Multiple sensors would require 
additional beamforming to factor in the arrival phase differences that would vary 
according to the sensor separations. 
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c. 


MEASUREMENT OF TRANSFER FUNCTIONS 


1. Setup 

The set up was performed in the anechoic chamber to minimize 
interference from reflections. The Microflown USP was rigged up in the array 1 
holder, and all three sensors were capped with their respective foam windshields, 
as they would be during data collection. 

2. Procedure 

The transfer functions for each of the three velocity sensors (BLUE, 
GREEN, RED) had to be determined separately. Therefore, the sensor had to 
be rotated each time so that the MRA of the each velocity sensor faced the 
source speaker. Broadband white noise was generated using a signal generator. 
Data from all sensors were collected via the National Instruments (Nl) cDAQ- 
9172 and processed in Matlab. 

3. Computation 

The required transfer function estimates are then generated in the 
frequency domain by taking the ratio of the cross-spectral density estimate 
between pressure (channel one) and the particle velocity channel (channel n) 
and the auto-spectral density estimate of the particle velocity channel: 



(x,(k)x !(*)) 

(*„(*) Km) 


n = 2 , 3,4 


(2.9) 


where H n (k )is the transfer function of k ,h frequency bin, X n {k ) is respective 
digital Fourier transform of the time-domain signals. 


1 This paper’s work only uses data from a single Microflown USP, as is sufficient to 
determine DOA as previously mentioned. However for practical reasons as well as to facilitate 
future work, the hybrid array holder made by Caulk (2009) is retained in the calibration process. 
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Figure 2 shows the plot of the transfer functions obtained for Sensor 323. The 
coherence subplot indicates that the usable frequencies begin from 
approximately 200 Hz onwards. 


Coherence - Sensor 323 



Phase Difference Sensor 323 



Frequency (Hz) 


Figure 2. Plot of transfer functions for Sensor 323 


7 



















THIS PAGE INTENTIONALLY LEFT BLANK 


8 



III. EQUIPMENT AND SETUP 


A. HARDWARE 

Appendix A provides a list of the hardware used. Table 1 lists some 
important configurations used for this thesis. 


Throughout entire study 

Microflown Four Channel Signal 

Conditioner 

- Gain set to‘HIGH’ 

- Correction mode set to ‘OFF’ 

Deriving Transfer Functions 

AUSTIN AU-15G Amplifier 

- Overdrive set to ‘ON’ 


- Volume set to ‘MAX’ 

HP Signal Generator 

- White Noise selected 


- Vpp set to ‘5V’ 

Windshield Testing 

Kestrel 4000 Handheld Anemometer 

- Data recording mode at 2-sec 

interval 


Table 1. Hardware configurations 


B. SOFTWARE 

Previous work done by Caulk (2009) had sensor data recorded via 
LabView, and subsequently post-processed in MatLab. This precluded real-time 
feedback for in-situ evaluation and adjustments if necessary. To overcome this, 
the latest Data Acquisition Toolbox within Matlab (version 2010b and above) was 
utilized to directly acquire data from the cDAQ. 


9 












The detailed component requirements are: 

• Matlab version 2010b and later 

• Data Acquisition Toolbox 

• Nl DAQmx 9.4 ( http://ioule.ni.eom/nidu/cds/view/p/id/2604/lanq/en ) 

This method also offered additional benefits such as programmatically 
setting the sampling frequency, auto-detection and processing of signals above 
threshold strength, etc. The downside however is that the responsiveness is 
dependent on the computer power available. This study used a Dell Precision 
M4500 laptop (Core i7 1,60GHz, 4 Gb RAM, 64-bit Windows 7) and was able to 
run the code continuously at a refresh rate of 0.5s or less. 

C. EXPERIMENTAL SETUP 

1. Anechoic Chamber 

For preliminary tests of the algorithms in the anechoic chamber, all of the 
equipment used for deriving the transfer functions were retained, except for the 
source. A bookshelf speaker was used, driven by a signal generator at 
monotones of various frequencies. Tests were conducted at varying locations 
respective to the sensor, as well as at different frequencies within the 2 kHz- 
Nyquist limit. 

A second set of tests also included ‘clapping hands’ and ‘humming’ while 
walking around the sensor array. Attempts were also made to test the system 
using mini remote-controlled helicopters, although it was observed that they did 
not generate enough signal strength. 

2. Field Measurement 

The primary purpose of this study was to acoustically detect and track 
aircraft. The chosen location was therefore in proximity (~1 mile to nearest edge 
of runway) to nearby Monterey Airport and positioned such that arriving and 

departing aircraft would almost fly overhead for maximum signal strengths. The 
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apparatus was mounted on a stand and placed on a concrete path next to a 
grass patch, approximately 30 meters from the nearest building (Hermann Hall, 
Naval Postgraduate School). The foam windshields were used on all sensors to 
minimize wind noise. Figure 3 shows a bird’s eye view of the field. 



Figure 3. Field location with respect to Monterey Airport’s runway 
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IV. SIGNAL PROCESSING 


This section discusses some of the programming and signal processing 
considerations. Key concepts and methodologies will be highlighted, together 
with notes about what worked and what did not. 

A. SAMPLING FREQUENCY 

The National Instruments Data Acquisition (DAQ) equipment uses a 
master internal timebase of 13.1072MHz. The sampling rate is controlled by the 
first installed module with a maximum of 51.2kHz. All other sampling rates must 
be integer divisions of the maximum according to the formula: 

/< = ^/256 ( 4 . 1 ) 

n 

where f M is the master clock rate and n is an integer from 1 to 31 (National 
Instruments, n.d). 

For this thesis, a sampling frequency of 4267Hz was chosen (using n = 12 
and rounded off) to achieve a Nyquist frequency of about 2 kHz. This was 
adequate for the application to commercial aircraft with noise signals typically 
ranging from 200 to 500 Hz. 

B. MATLAB / DAQ-MX CONFIGURATION 

Matlab’s Data Acquisition Toolbox works with Nl DAQmx to import data 
from the Microflown USP, via the Nl cDAQ in real time. The process is generally 
straightforward and simple to use. The following discusses some salient points 
during set up: 
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1. Listening Mode 

Two modes are available; Foreground or Background. The latter was 
chosen since it allowed simultaneous operations to be carried out while the data 
was streaming. Foreground mode presumably would be more suited to 
applications requiring higher degrees of control, such that data processing 
commences only upon the completion of data capture. 

2. Data Stream / Buffer 

It was observed that the data stream in Background mode enters a 
memory buffer, and when deemed ‘available’ the programmed function will be 
invoked. It appeared, however, that the ‘Data Available’ criterion was not 
consistent either by data size (number of data points) or time. A preliminary 
conclusion was that this refresh rate was a function of overall computing speed 
and data capture parameters, possibly including the prescribed sampling rate 
and even the USB transfer rate. The variance of each sequence was however 
not more than 10% each time. For example, in one instance the buffer could 
provide 590 data points, the next 605, then 610. 

Nonetheless, the subsequent data preparation and beamforming expects 
a consistent and pre-determined data length. As a workaround, the data stream 
was cumulatively assigned to a growing vector, and the latest N length of data 
would be used each time. The length N would be a function of the desired 
resolution (longer means better frequency resolution after FFT), noise reduction 
potential (more data allows noise reduction through averaging) as well as target 
motion resolution (shorter means a tighter snap-shot in time). This thesis 
adopted the following: 

Length of data processed each time = 2 10 = 1024 

Given a sampling frequency, f s of 4265 Hz, 

f 

Frequency resolution = ^- + 1= 4.16 Hz 

N 
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Time duration 


1024 

~T 


= 0.24 s 


3. Real-time processing and feedback 

This is a key component of this thesis. To allow at least near real-time 
feedback, the data stream needs to be processed and its results presented 
almost immediately. The time allowed for this would naturally vary across 
applications, although in general the faster change in target position (as 
observed by sensor), the faster refresh rate necessary. 

Depending on the available computer power, no additional setup may be 
necessary. The Dell M4500 laptop, however, was not able to process every data 
stream quickly enough. The alternative was therefore to only process data 
streams at intervals, say every x instances of ‘Data Available’. In this case, 
processing only every 4th data stream resulted in a sustainable refresh rate of 
less than 0.5 seconds. Each time only the latest N bits of data were used to 
ensure the most up-to-date data was being processed. 

4. Saving Data Streams 

Notwithstanding the above, the integrity of the entire input data stream is 
preserved in the operating environment, and is saved to file for further processing 
if needed. This allows detailed analysis which may not have been possible 
during collection. 

C. FREQUENCY DOMAIN ANALYSIS 

Once the data is properly extracted, most of the processing occurs in the 
frequency domain. This allows the beamforming algorithms to target specific 
frequency or frequency bands of interest. The following discusses some 
techniques employed to enhance results: 
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1. Bandpass Filtering 

The power spectrum of the incoming data can be evaluated in real-time, 
and a suitable pass band can be determined through a combination of prior 
programming bias, and/or rules such as maximum energy detection. This thesis 
employed a pre-determined pass band of 300-600 Hz, which corresponded to 
most signals emitted by commercial aircraft. Since the data has already been 
transformed to the frequency domain, this is a simple process of discarding data 
outside the pass band. 

2. Peak Detection 

Within the pass band, a peak detection algorithm determines the peak 

frequencies for further evaluation. The algorithm compares each bin value to 
that of its neighbors, and classifies those above the preset threshold as peaks. 
Threshold settings are adjustable for further refinements. 

3. Beamforming 

This is the main process of taking the dot products between the MRAs of 

various particle velocity channels and the unit vector in the steered direction. As 
expected, this part is the single largest computational bottleneck of the entire 
program. The more frequencies selected for beamforming and the finer the 
angular resolution, the longer it takes and slower the overall refresh rate. 

Following the peak detection routine, various approaches are possible: 

a) Single maximum peak (plus buffer above and below) 

b) All detected peaks (plus buffer above and below) 

c) Any arbitrary pass band 

Initial configurations to achieve the less-than-0.5 second refresh rate 
adopted Option (a). Option (b) proved to be marginally slower, and Option (c) is 
highly intensive and more suited to be a post-processing approach. 
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4. Angular Resolution 

The code also allows the angular resolution to be adjustable, and the 
obvious tradeoff would be between that of accuracy and speed. This was set to 
be 1-degree increments for the field experiment, and 5-degree increment for the 
tests conducted in the anechoic chamber. 

D. VISUAL DISPLAY FOR FEEDBACK 

Last but not least, the results of beamforming has to be presented for 
feedback. Caulk’s work included intensive graphics such as 3D surface plots 
and full-color imagery. In order to achieve near real-time feedback, this had to 
be replaced by basic 2-D plots (c|) - 0), as well as a 3-D directional cone display. 
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V. RESULTS AND ANALYSIS 


A. IN ANECHOIC CHAMBER 

1. Single Source, Varying Frequency 

Driving a single speaker at various frequencies within the Nyquist limit 
gives clear and distinct results. Utilizing the single peak detection (Option (a)) 
also allowed for fast refresh rates and almost instantaneous feedback. The 
system performed equally well for most frequencies, with the lowest workable 
frequency found to be approximately 300 Hz. This is slightly higher than the 200 
Hz suggested by the earlier coherence plots. The following tests will be 
presented here: 

a) 1000 Hz along positive Y-axis (Figure 4) 

b) 1000 Hz along positive Z-axis (Figure 5) 

c) 300 Hz along positive Y-axis (Figure 6) 

The figures show the detected peak in frequency domains, the 2-D plot, 
as well as the 3-D cone plot intended as visual feedback in real-time. The figure 
titles indicate the detected frequencies, as well as source positions. 
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Figure 4. Single Source at 1000 Hz along positive Y-axis 




Figure 5. Single Source at 1000 Hz along positive Z-axis 
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St ogle-Sided Amplitude Spectrum of y(t> 



Figure 6. Single Source at 300 Hz along positive Y-axis 

While the 1000Hz tone showed excellent tracking along both Y and Z 
directions, inaccuracies started showing in the tracking of the 300Hz tone. When 
placed along the positive Y-direction, its detected location was at cj) = 85° and 0 = 
110° (5° and 20° error respectively). 

2. Two Source, Different Frequency 

Two speakers were placed at approximately the same height of and 
distance away from the sensor array, separated by approximately 50 degrees in 
azimuth (in the x-y plane). One speaker was driven by a 1000Hz monotone, the 
other a 300Hz. The system managed to sequentially localize each source, 
depending on the frequency selection (choice of peak). 
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Figure 7. Two sources (300 Hz and 1000 Hz), 300 Hz detected. 



Figure 8. Two sources (300 and 1000 Hz), 1000Hz detected. 

Figure 7 and Figure 8 show the correct detection of both monotones 
emitted simultaneously, albeit with some inaccuracies. The heights of both 
sources were the same as the previous test (i.e. 9 = 90°), but the 300 Hz tone 
was detected at 0 = 125°, while the 1000 Hz tone was detected at 0 = 105°. 
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While it was not done, the system could easily be modified to 
simultaneously locate and display both source locations transmitting at different 
frequencies. This would be done via plot-holding in Matlab, and overlaying the 
two distinct sources as detected by the peak-finding algorithm. 

3. Two Source, Same Frequency 

As expected and presented by Hochwald and Nehorai, the basic 
beamforming algorithm employed is unable to distinguish between two sources 
transmitting at the same frequency or at overlapping frequency bands. A method 
known as Multiple Source Classification (MUSIC) is widely known to be able to 
do that. Given n vector sensors, the theory hypothesizes the detection of 4n-2 
uncorrelated sources. A single vector sensor should therefore be capable of 
distinguishing two sources transmitting at overlapping frequencies. 

B. FIELD RESULTS 

1. Original Setup 

A variety of aircraft were tracked in this session, the two most distinctively 
different ones are the twinjets and the twin turboprops. 
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LOFARgram of Arriving E55P Twinjet 


Figure 9. LOFARgram of arriving twinjet 

Figure 9 shows the LOFARgram of an arrival twinjet. A few things can be 
seen immediately from the plot: 

a) Target SNR is high and clearly detectable above ambient noise 

b) Most ambient noise occur at <200 Hz 

c) Clear Doppler shift 

d) Distinct Lloyd’s Mirror pattern due to interference from reflections 
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LOFARgram of Arriving E55P Twinjet 



Frequency (Hz) 


Figure 10. LOFARgram of arriving jet with detected peak frequencies 

Figure 10 shows the same LOFARgram, as well as the detected peak 
frequencies by the peak-finding algorithm, which would be subsequently used for 
beamforming. Within the specified range of 300-600 Hz, the peak-finding 
algorithm worked well to identify the main signal frequency most of the time. The 
detected peaks also helped to clearly identify the Doppler shift. 
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Arriving Jet (Real-time) 
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Figure 11. Real-time cj) - 0 plot for arriving jet 

The real-time tracking, presented as a 2-D cj) - 0 plot is shown in Figure 11 
and was clearly not an accurate representation of the actual flight path. Since 
the plane was landing, it can be assumed that it was on a straight course. Even 
with a changing altitude, the actual cj) - 0 plot should be relatively smooth, with 
small incremental changes in the polar and azimuthal angle. 

Since the full data was recorded, post-analysis was possible. The range 
of beamforming options as earlier presented was tested, and the best result was 
produced by beamforming across a wider range of frequencies. The following 
series of plots (Figure 12) illustrate the incremental improvements. 
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Arriving Jet (Real-ome) Arriving Jet (300-600Hz) 




Amving jet (300-1000 Hz) Amving Jet (300 -1500 Hz) 




Figure 12. Significant improvements in tracking by beamforming across wider 

range of frequencies—Arriving Jet 


Average beamforming times per data point/cycle are appended 
accordingly. Beamforming across more frequencies essentially perform 
averaging, and reduces the effect of noise and other interferences. In particular, 
a vector sensor is most sensitive to reflections that alter the amplitude and phase 
information of the various channels, and averaging may have mitigated that. 

Instead of a selected few frequency bins, the beamformer now computes 
many more and then sums them cumulatively, before presenting the global 
maxima for cfi and 0 respectively. This increased computational load would have 
to be balanced with the desire for system responsiveness. 
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B. DEPARTING TURBOPROP 


LOFARgram of El20 Turboprop 



Frequency (Hz) 


Figure 13. LOFARgram of departing turboprop 

Figure 13 shows the LOFARgram of a departing turboprop. As compared 
to the jet, the broadband signal strength is weaker but still clearly above ambient 
noise levels. There is still a detectable, albeit smaller Doppler shift due to the 
apparently lower landing speed. A Lloyd’s Mirror pattern can still be seen, but is 
less obvious due to the smaller signal amplitudes. This is consistent with the in- 
situ observation that the jets sound louder than the turboprops. 
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Figure 14. LOFARgram of departing turboprop with detected peak frequencies 


Figure 14 shows that the peak-finding algorithm is now struggling to ‘lock’ 
onto a consistent source. Within the same specified range of 300-600 Hz, the 
presence of similar-amplitude harmonics appear to be challenging for a single 
source to be tracked closely. It can also be seen that some high-energy peaks 
are occurring at around 200 Hz. While expanding the peak-finding algorithm 
would produce a better ‘lock’ on these emissions, earlier tests had shown the 
beamformer to only perform decently at frequencies above 300 Hz. 

That said, all the detected peaks are without doubt transmissions 
originating from the target, and not elsewhere. Theoretically, a vector sensor 
should therefore still be able to resolve the target position from these frequency 
peaks. 
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Figure 15. Real-time (j) - 9 plot for arriving turboprop 


As seen in Figure 15, the tracking of the turboprop in real-time was also 
not accurate. Various other beamforming options, including the proven method 
of summing across more frequencies in the arriving jet did not produce similar 
improvements. This is shown in the following series of plots (Figure 16). 
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Departing Turboprop (Real-time) Departing Turboprop (300 - 600 Hz) 




Departing Turboprop (300 - 1000 Hz) Departing Turboprop (300 • 1500 Hz) 




Figure 16. Negligible improvements by beamforming across wider range of 

frequencies—Arriving turboprop 


A few reasons could explain this, the most likely of which is that the 
highest energy emissions occurred below the minimum working frequency of the 
Microflown USP. The higher harmonics may not be adequately strong to provide 
the correct beamforming, or were severely affected by reflections. 

2. Revised Setup—with Foam Box 

Following the results from the first field setup, an attempt was made to 
minimize bottom reflections. A cardboard box, lined with sound-absorbent foam 
was used to shield the bottom half of the sensor. To simplify matters, the sensor 
array was also replaced by a single Microflown USP on an equipment stand. 
Figure 17 shows a picture of the revised setup. 
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Figure 17. Revised setup of single Microflown USP and anti-reflection foam 

box still need picture 

Since the transfer functions had been derived for the sensor array, the 
new setup technically had to be re-calibrated by performing the same procedure 
in the anechoic chamber. A quick check in the chamber however verified that the 
old transfer functions were still valid, and the new setup could still track sources 
accurately. The re-calibration was therefore not performed in the interest of time. 
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Arriving Twinjet 


LOFARgram ot Arriving Twinjet (with foam box) 
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Figure 18. LOFARgram of arriving twinjet (new setup with foam box) 
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Figure 19. LOFARgram of arriving twinjet (original setup) 

Figure 18 shows a markedly different signal profile of a twinjet, with 
a foam box now under the sensor. Compared to the jet profile using the original 
setup, shown in Figure 19, the Lloyd’s Mirror pattern has been reduced 
substantially, though still noticeable. Other aspects such as the Doppler shift and 
overall signal amplitude also seem to have been attenuated to some degree. 
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Figure 20. Real-time cj) - 0 plot for arriving twinjet (new setup with foam box) 


Figure 20 shows that despite the modifications, the real-time 
tracking is still not accurate. Post-analysis however show greater incremental 
improvements when summing across more frequencies. 
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Arriving jet with foam box (realtime) Arriving jet with foam box (300 - 600 Hz) 




Arriving jet vMth foam box (300 -1000 Hz) Arriving jet v*th foam box (300 • 1500 Hz) 




Figure 21. Faster improvements by beamforming across wider range of 
frequencies—Arriving Jet (with foam box) 


Figure 21 shows that the tracking improved substantially faster with 
beamforming across more frequencies, as compared to the previous setup. The 
tracking however, while relatively linear and smooth, did not depict the actual 
flight path. With 0 reaching 160° at the closet point of approach (t() = 0°), this 
tracking meant that the flight path was only at a 20-degree elevation from the 
sensor, when it clearly flew almost overhead. Closer analysis revealed that the 
the vertical particle velocity sensor (GREEN, x-axis) was recording a less-than- 
expected signal amplitude. This produced the erroneous tracking of a low-flying 
target. 

It is unclear what may have caused this, but suspicions point to the 
newly added foam-lined box underneath the sensor. The presence of a faint 
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Lloyd’s Mirror pattern shows the box is not completely opaque to the signals, it 
may have in other ways interfered and contributed to the specific attenuation to 
the GREEN channel. Preliminary research showed that the wedge foam used 
was not particularly effective in absorbing the low frequencies of 250-500 Hz. If 
the GREEN signals were artificially boosted by a applying a scalar multiple, the 
resultant tracking approaches the true flightpath, as shown in Figure 22. 


Arriving jet with foam box (300 • 600 Hz) 



Arriving jet with foam box - GREEN channel x 5 



Amvirig jet with foam box - GREEN channel X 10 



Arriving jet with foam box - GREEN channel x 20 



Figure 22. Improved tracking by boosting GREEN signal by 5, 10 and 20 times 

(beamformed300-600 Hz) 
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b. Arriving Turboprop 



Figure 23. LOFARgram of arriving Turboprop (new setup with foam box) 




Figure 24. LOFARgram of departing Turboprop (original setup) 


Figure 23 shows the signal profile of an arriving turboprop with the 
foam box. Compared to that using the original setup in Figure 24, the Lloyd’s 
mirror pattern is noticeably suppressed, especially towards the higher 
frequencies. Signal strength before and after CPA also appeared to have been 
reduced to some degree, resulting in shorter frequency ‘streaks’. 
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Arriving turboprop (Real-time) Amving turboprop (300 - 600 Hz) 




Amving turboprop (300 -1000 Hz) Amvtng turboprop (300 -1500 Hz) 




Figure 25. Real-time t)> - 0 plot for arriving turboprop (new setup with foam box) 

Figure 25 shows that the real-time tracking again is not accurate. 
Beamforming across more frequencies also did not give much improvement. 
Compared to those using the original setup in Figure 16, the tracks, however, 
were smoother and bore greater resemblance to a realistic flight path. 

The 300-1500 Hz plot in particular was interesting. This implied 
that the plane flew past on the other side of the sensor, away from the buildings 
(0 < 90°), and at an elevation of 20°. This was previously unseen, and could be 
due to the reflections off the wall of the cardboard box. Furthermore, the 
modifications would still not address the fact that some high-energy emissions 
from the turboprops are occurring below 300 Hz. These would continue to 
present problems for the Microflown USP since we were unable to produce loud 
enough low frequency noise to establish good transfer functions that far down in 
frequency. 
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3. Analysis of Lloyd’s Mirror Pattern 
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Figure 26. LOFARgram of a arriving jet with distinct Lloyd’s mirror pattern 


Despite reflections being undesirable for vector sensors in general, 
additional useful information can be derived from careful analysis. This section 
explains. The salient features, as indicated on Figure 26 are: 

• Doppler shift ‘streaks’ 

• Locations of maxima at CPA 

• Gradients of frequency peaks 
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a. Target Velocity 

This is the direct application of the Doppler shift equation valid 
where the contact speed is small compared to the speed of sound: 

f - f 

Equation Section 5 v = c— —— (5.1) 

fu + fl 

where v is the velocity of target, c is the speed sound in air (343 m/s) 
and f u , / ; are the upper and lower frequencies of the Doppler streak. 


LOFARgram of arriving jet (Lloyds's mirror pattern) 
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Figure 27. Deriving target velocity from Doppler shift 

Choosing the highest frequency visible Doppler streak in Figure 27 
for the lowest minimum errors, f u = 1770 Hz, and f u - 1204 Hz. 


v = 343 


1770-1204 
1770 + 1204 


= 65 ms '(12 6 knots) 
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b. Target Altitude 

Assume the sensor has a height d , and the aircraft is flying in the 
y-direction at speed v and altitude h . The horizontal range along the x-axis at 
CPA is i? 0 and occurs at time t - 0, as illustrated in Figure 28. 



Figure 28. Schematic of flypast 

The exact solution for the path length difference between the direct 
and reflected paths is: 


R reflected ~ R ,direct ~ + ^) + ^0 + " >/(^ “ <0 + ^0 + 


2.2 

V t 


(5.2) 
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As /-» oo, v 2 f 2 becomes large compared to [h+df+R^ and the 
path length difference can be approximated as: 


A r ■ 




2 

2L2 


A r = vt 


l4 1 (h + d) +R~ j 4 l jh-d) +Rp ^ 


O 2.2 

2 v t 


o 2.2 

2 v t 


(5.3) 


Expanding the terms and simplifying results in: 

2 hd 


Ar = ■ 


vt 


(5.4) 


Assuming a rigid boundary condition for the reflection, when the 
path length equals an integral number of wavelengths, you get maxima in the 
LOFARgram. This condition is given by: 


A r = 


2 hd 
vt 



Solving for time and differentiating with respect to frequency gives: 

dt ~ 2 hd 
df n ncv 


h 


dt (ncv ' 

df n l 2d j 


(5.5) 


(5.6) 


So at times far from CPA, the time-frequency plot of the maxima 
are linear with a slope which decreases by integer multiples from the steepest 
slope for n = 1. Measurement of the slope would therefore yield the altitude of 
the aircraft. 
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LOFARgram of arming jet (Ltoyds's mirror paaem) 
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Figure 29. Deriving target altitude from gradient of frequency peaks 


Applying this result to the example above, we first determine the n 
= 1 peak at approximately 182 Hz (Figure 29). This can be checked by verifying 
the subsequent maxima are indeed integer multiples of this frequency. For 
example, 2nd maximum is 365 Hz (n = 2.00), 3rd maximum is 552 Hz (n = 3.03). 

Next we choose a gradient of frequency peaks. Recall the use of 
t -» oo in simplification earlier. It is therefore best to choose the visible frequency 
gradient that is furthest away from CPA. 

In this example, the n = 3 maxima was used. Using d = 0.9m and 

n = 1, 


h = 


6.971-5.056 ^ (3)(343)(65.2784) ' 
1091-733.8 [ 2(0.9) J 


= 200 m 
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c. 


Horizontal Distance at CPA 


At CPA, t = 0. (5.2) then simplifies to: 


Ar cPA=\j( h + d ) + R q- (h-d ) +R { 


This path length difference produces maxima at: 


(5.7) 


f n =~ = 7 (5-8) 

yj(h + d) +R 2 - - x j(h -d) + R 2 

Using the previously estimated altitude/?, and the frequencies of 
the maxima observed at CPA, the exact solution for /? 0 can be found using a 
standard root finding algorithm. If the horizontal range at CPA is small compared 
to the altitude (R g « h-d), as in this case of overflying targets, R 0 can also be 
estimated explicitly. 


A r = ^(h+df + R 2 0 -yj(h-df+I% 


Ar = (/i + d) + " 


(, h+df 


(h-df 


Assuming R « h-d), this expression can be simplified to: 


Ar = (h + d) 


f 1 n2 A 

1 + 1 ^ 


V 


(h+dy 


-(h-d) 


f , % \ 

t , 1 * 


V 


(h-d)- 


(5.9) 
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Expanding the terms and simplifying gives: 


A r = 2 d + 


1 *p- 

2 (h + d) 


1 Rq 

2 (h-d) 


= 2 d + 


1 

R o 

1 

R l 

2 , i 

( d^ 

i 2,i 

r d\ 

h \ 


I h\ 



l h) 


l h) 


(5.10) 


Now assuming h»d (200m vs 0.9m in this case), we can further simplify this 
into: 


, 1 /< 

f 

d "j 

' 1 

( d 'l 

R 2 

f R n) 

Ar = 2d+ -5. 

1- 




= 2d - ird = d\ 

2 - f 

2 h 1 

V 

Tr) 

1 2 A 1 

l h 

1 h 2 

l h 2 ) 


(5.11) 


From (5.5), 


a nc J 

A r = — = d 


fn 


2 _ R 1 

h 2 


R=h 




(5.12) 


Applying to our example, and using the first maximum at 182 Hz, 
h = 0.95 m (using 0.9 m would result in an imaginary root), d = 200 m, 


R = 200 2 — 


(1)(343) 

(182)(0.95) 


= 25 m 


The exact solution using (5.8) and Matlab ‘roots’ function gave an answer of 25.6 
m, thus validating the prior assumptions. 
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VI. CONCLUSION 


The use of a single acoustic vector sensor in tracking commercial aircraft 
in final landing approaches has been successfully demonstrated. The 
advantages of a relatively small package with wideband capability could be 
substantially exploited in mobile systems. 

The specific setup used however, did not manage to do so in real-time as 
desired, which is essentially a function of computing power. Algorithms that are 
more efficient could potentially improve performance. 

The results revealed that vector sensors, since relying heavily on the 
amplitude and phase information between the particle velocity and pressure 
channels, are especially sensitive to reflections. The contrast between in¬ 
chamber and field testing illustrate the adverse effects of bottom reflections. 
Curbing such interferences would therefore be an utmost priority in the 
employment of vector sensors in general. Nonetheless, it was also illustrated 
that reflections can become a useful source of information when properly 
analyzed. 

In seeking an eventual capability of real-time tracking, continued efforts 
could be made in various directions: 

• Tracking of higher-frequency (500 Hz and above) targets which operate 

well clear of the minimum working frequency 

• Using acoustic-absorbent material designed for frequency band of interest 

• Re-positioning sensor to minimize bottom reflections 
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APPENDIX A. LIST OF EQUIPMENT 


Microflown Holland Four Channel Ultimate Sound Probe 
(USP)—Model UT0901-23/24 (Sensor 323/324) 



Microflown Holland Four Channel Signal Conditioner— 




Model E0901-23/24 


National Instruments CompactDAQ USB Chassis 
Model Nl 9172 


National Instruments Sound and Vibration DAQ Module- 
Model Nl 
9234 (3 total) 


AUSTIN AU-15G Amplifier 


Hewlett Packard Function Generator—Model 33120A 


Kestrel 4000 Handheld Anemometer 


Dell Precision M4500 Laptop Computer 
(Core i7 1.60 GHz, 4 Gb RAM, 64-bit Windows 7) 


- **** • 
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APPENDIX B. COMPUTER CODE 


A. MAIN SCRIPT 

main.m 

%This is the main script to run in real-time, and saves the entire data 
%for subsequent processing and analysis. User to specify: 

%Fs - sampling frequency in Hz 

%duration - collection duration in seconds 

clear all 
close all 
clc 

%Specify parameters 

Fs = 4265; %samplng frequency 

duration = 5; %seconds 

%use acquireData 4ch or acquireData_12ch as appropriate 
data = acquireData_4ch(Fs, duration); 

%save data to file with system timestamp 
dt = datestr(now, 'mmmdd HHMMSS' ); 
ft = strcat(' .\Nfloor\' ,dt); 
save (ft, 'data' ) 
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B. KEY FUNCTIONS 


acquireData 4ch.m 


function [data] = acquireData_4ch (Fs, duration) 

%Starts the data acquisition process using Matlab's Data Acquisition 
%Toolbox, and via NI DAQmx drivers. 

%Inputs from calling function: 

%Fs - sampling frequency in Hz 

%duration - acquisition duration in seconds 

%Variables need to be global so to be used across various functions 
global data 
global Fs 

%Create generic session for NI devices 
s = daq.createSession ( ' ni ' ); 

%Setting parameters using provided inputs 
s.Rate = Fs; 

s.DurationlnSeconds = duration; 

%Add device and channel/s. 

%If unsure, type 'daq.getDevices()' to see what's connected and 
available 

s.addAnaloglnputChannel( 'Devi ',0:3, 'voltage' ); 

%The previous steps take a few seconds, and may affect timely 
recording. 

%The following pause inserts a break, so recording starts almost 
%immediately after a key stroke here, 
pause () 

%Start background listening, and specify function to process data when 
%available 

Ih = s.addlistener( 'DataAvailable' ,@plotData); 
s.startBackground(); 
s.wait() 
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acquireData 12ch.m 


function [data] = acquireData_12ch (Fs, duration) 

%Starts the data acquisition process using Matlab's Data Acquisition 
%Toolbox, and via NI DAQmx drivers. 

%Inputs from calling function: 

%Fs - sampling frequency in Hz 

%duration - acquisition duration in seconds 

%Variables need to be global so to be used across various functions 
global data 
global Fs 

%Create generic session for NI devices 
s = daq.createSession( 'ni' ) ; 

%Setting parameters using provided inputs 
s.Rate = Fs; 

s.DurationlnSeconds = duration; 

%Add device and channel/s. 

%If unsure, type 'daq.getDevices()' to see what's connected and 
available 

s.addAnaloglnputChannel( 'cDAQlModl ',0:3, 'voltage' ); 
s.addAnaloglnputChannel( 'cDAQlMod2 ',0:3, 'voltage' ); 
s.addAnaloglnputChannel( 'cDAQlMod3 ',0:3, 'voltage' ) 

%The previous steps take a few seconds, and may affect timely 
recording. 

%The following pause inserts a break, so recording starts almost 
%immediately after a key stroke here, 
pause () 

%Start background listening, and specify function to process data when 
%available 

Ih = s.addlistener( 'DataAvailable' ,@plotData); 
s.startBackground(); 
s.wait() 
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plotData.m 


function plotData(src,event) 

Called by 'acquireData_4ch.m' or 'acquireData_12ch.m' when data is 
available. Appends entire data stream to global variable 'data' 

persistent tempData; 
persistent i; 

global data 
global Fs 

if (isempty(tempData)) 
tempData = []; 


end 

data = [data;event.Data]; 


%Specify data chunk length for subsequent processing 
N = 2 A 10 ; 

NFFT = N/2; 

i = floor(length(data)/N); 


j = 0; 

if i > 0 S i > j 

start = (i—1)*N+1 
finish = i*N 

y = data(start:finish,:); 

%Process every 'x' available data chunks to adjust system 

response 

%time. Balance between refresh rate and system lag. 
x = 4; 

if (mod(i,x)==0) 

beamform 4ch(Fs,N,NFFT,y); 

end 


%Plots single-sided amplitude spectrum for feedback in real¬ 
time 

L = length(y); 

NFFT2 = 2 A nextpow2(L); 

Y = fft(y,NFFT2)/L; 
f = Fs/2*linspace(0,1,NFFT2/2+1); 

figure (1) 

y2 = 2*abs(Y(1:NFFT2/2+l)); 
plot(f, 2*abs(Y(1:NFFT2/2 + l))) 
axis([0 2000 0 50*mean(y2)]); 
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time = num2str(i); 

time = strcat({ ' '}, time, {' sec'}); 

figtitle = strcat (' Single-Sided Amplitude Spectrum of y(t)', {' 

' }, time); 

title(figtitle) 

xlabel (' Frequency (Hz) ') 

ylabel ( ' |Y(f) | ' ) 


j = i 

end 

end 
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beamforme 4ch.m 


function beamform 4ch (Fs,binN,N,bindata) 

%use tic & toe to time beamformer cycle 
tic 

global Fpoints targetPosition targetphi targettheta time_taken; 


%DAQ channel assignments 

prsl = 1; %sensor 324 pressure 

blul = 2; %sensor 324 blu velocity Z axis 

grnl = 3; %sensor 324 grn velocity X axis 

redl = 4; %sensor 324 red velocity Y axis 

thetainc = 1; %angle increment degrees 
phiinc = 1; %angle increment degrees 
numsens = 1; %total number of sensors 
numel = 4; %number of elements per sensor 


%calculated parameters and constants 
k = linspace(0,N-l,N); %baseline FFT bin numbers 
d = 0.172; %array element spacing meters - measured 172mm ± 2mm 
c = 340; %speed of sound meters/sec @ 20°C sea level 
maxchan = numsens * numel; %total channels 
x = 1; %X axis designator 
y = 2; %Y axis designator 
z = 3; %Z axis designator 

%Vector Sensor Angles relative to Array Apparatus - radians 

thetabl = 0; %323 blul 0° 

phibl = pi/2; %323 blul 90° 

thetagl = pi/2; %323 grnl 90° 

phigl = 0; %323 grnl 0° 

thetarl = pi/2; %323 redl 90° 

phirl = pi/2; %323 redl 90° 


sSensor Element Unit Vectors 


ue(prsl,:) =[00 0]; %323 prsl dummy 
[omni] 

ue(blul,:) 
cos(thetabl)]; 

ue(grnl,:) 
cos(thetagl)]; 

ue(redl,:) 
cos(thetarl)]; 


= [sin(thetabl)*cos(phibl) 
%323 blul [0 0 1] 

= [sin(thetagl)*cos(phigl) 
%323 grnl [1 0 0] 

= [sin(thetarl)*cos(phirl) 
%323 redl [0 1 0] 


[000] - change in udot to 
sin (thetabl)*sin(phibl) 
sin(thetagl)*sin(phigl) 
sin(thetarl)*sin(phirl) 


%Steered Unit Vectors 

thetas = 0:deg2rad(thetainc):pi; %polar steer angle array- radians 
phis = -pi:deg2rad(phiinc):pi; %azimuthal steer angle array upper 
hemisphere - radians 

%phis = -pi/2:deg2rad(phiinc):pi/2; %azimuthal steer angle array 
upper hemisphere - radians 
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m 


= length(thetas); %total number of angles in theta direction 
n = length(phis); %total number of angles in psi direction 
p = length(k); %vector length in k direction 

u = sin(thetas') * cos(phis); %steered direction cosine x - m x n 
matrix 

v = sin(thetas') * sin(phis); %steereddirection cosine y - m x n 
matrix 

w = cos(thetas') * ones(l,n); %steereddirection cosine z - m x n 
matrix 

%element unit vector & steer angle dot product array 
%calculate one udot per element, same for all sensors, 

%use first four channels of ue matrix 

%vectorize dot products for faster multiplication in beamformer 
udot = ue(prsl:redl,x)*u(:)' + ue(prsl:redl,y)*v(:)' + 
ue(prsl:redl,z)*w(:)'; %l-4 x m*n 

%omni-directional pressure sensors - vector is "1" in all 
orientations and aspects 

udot(l,:) = ones(l,m*n); %1 x m*n 


Xc = zeros(maxchan,N); 

%hanning window to reduce sidelobes of signal 
H = hann(binN) * ones(1,maxchan); 
xwin = bindata .* H; 

%clear bindata H 

%convert binary signal data to frequency domain 
Xw = fft(xwin); %by columns binN x maxchan 

%truncate length of signal to match transfer function length N x 
maxchan 

X = Xw(1:N,:); 

X = X.'; %change to row data maxchan x N, non conjugate transpose 
clear xwin Xw %save memory 


load '.\txfunction 323.mat' 

Hplb_BLUE = Hp1B(1:N, : ); 

Hplb_GREEN = HplG(l:N,:); 

Hplb_RED = Hp1R(1:N, :); 

%Apply transfer functions 

Xc(1,:) = X(1,:); 

Xc(2,:) = X(2,:) .* Hplb BLUE'; %Hbl=Hpl*Hbl - pre multiply vel and 

prs Xfer funct 

Xc(3,:) = X(3,:) .* Hplb GREEN'; %Hgl=Hpl*Hgl - pre multiply vel 

and prs Xfer funct 

Xc(4,:) = X(4,:) .* Hplb RED'; %Hrl=Hpl*Hrl - pre multiply vel and 

prs Xfer funct 


%apply band-pass filter 
length(Xc) 
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startband = 300; 
stopband = 600; 

startbin = floor(startband*2*N/Fs) ; 
stopbin = floor(stopband*2*N/Fs); 

Xc2 = zeros(raaxchan,N) ; 

Xc2(:,startbin:stopbin) = Xc(:,startbin:stopbin); 

%Use peakfinder.m to determine source frequency/s 
[peakLocl, peakMag]=peakfinder(Xc2(1, :) ) ; 

peakLoc = peakLocl.*Fs./(N*2); %convert from bin numbers to actual 
frequencies 

%choose maximum peak within earlier specified passband 
[C,I] = max(peakMag); 

Fo = peakLoc(I) ; 

Fpoints = [Fpoints;Fo]; %store into global variable for appending 
onto LOFARgram 


%+- bins around design freq to truncate - processing bandwidth, 
ignore remainder 

Fplot = floor(Fo*(N*2)/Fs) + 1; %k bin for freq Fo 

HBW = Fo/100; % + - freq to offset signal from Fo - broad for live 
signals 

halfbandwidth = floor(HBW*(N*2)/Fs) + 1; %+- number of bins to 
offset signal 

kshort = k(:,Fplot-halfbandwidth:Fplot+halfbandwidth); %get 
relevant k values 

k = kshort; %transfer relevant k values to new k vector 

clear kshort 

p = length(k); %new vector length in k direction 


%select appropriate bins for beamforming 

%Option 1: Use detected peak frequency Fo +- buffer 

Xcshort = Xc(:,Fplot-halfbandwidth:Fplot+halfbandwidth); 

%Option 2: Use continuous range of frequencies for noise reduction 
average 

%Xcshort = Xc(:,startbin:stopbin); 

Xc = Xcshort; 

clear Xcshort %save memory 
Xc = conj(Xc); 


%%initialize beamformer%% 

%%S is sum of all elements and sensors - the beamformer%% 

S = zeros(m*n,size(Xc,2)); %final beam former output - m*n x p 
array 


%Beamformer output 
for q = 1:4; 

%pull element matrix from udot array 
%Theta m * Phi n vector 

udotc = squeeze(udot((mod(q+3,4)+1),:,:)); %m*n x 1 
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end 


sumcalc = udotc(:) * Xc(q,:); %m*n x p 
S = S + sumcalc; 


S = sum(abs(S),2); %sum along k "amplitude", collapse to single 
value in k direction 

S = reshape(S, [m n]); %mxn thetaxphi 2D matrix 

%remove SdB points outside 90% max value for clearer phi-theta 
plots 

indices = S > max(max(S)) .*0.9; 

S = indices.*S; 

%convert to dB 

SdB = 20*logl0(abs (S)); 

%Find target phi and theta by searching for maximum values 
[SdBmaxl,11] = max(max(SdB, [ ],1)); 

[SdBmax2,I2] = max(max(SdB, [ ],2)); 
phisdeg = rad2deg(phis); 
thetasdeg = rad2deg(thetas); 
targetphi = phisdeg(II); 
targettheta = thetasdeg(12); 
targetphirad = phis(II); 
targetthetarad = thetas(12); 

%Append onto global vector for final tracking plot 
targetPosition = [targetPosition; [targetphi targettheta]]; 

%Convert into coordinate system of 'Cone' subfunction 
x = sin(targetthetarad) * sin(targetphirad) * 20; 
y = cos (targetthetarad) * 20; 

z = sin(targetthetarad) * cos (targetphirad)* 20; 

Cone([0 0 0],[x y z],[0 l],10,'r',0,0,Fo,targettheta,targetphi); 


%Overhead Theta by Phi amplitude plot 
figure (5) 

%clims = [min(min(SdB)) max(max(SdB))]; 

%imagesc (rad2deg (phis) , rad2deg (thetas) , SdB, dims) ; 
imagesc(rad2deg(phis),rad2deg(thetas) , SdB) ; 

location = strcat({ ' \phi='}, num2str(targetphi), {' \theta='}, 

num2str(targettheta) ) ; 

figtitle = strcat ({' Detected source of '}, num2str(Fo, 1 %4.Of' ), {' 

Hz'], location); 

title(figtitle) 
axis([-180 180 0 180]) 
xlabel( '\phi (deg)' ) 
ylabel( '\theta (deg)' ) 
shading interp 

%{ 
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%3D rendering of sensor beam pattern. This allows visual 
verification 

of vector sensor's performance, which is has a maximum in the 
direction 

of source, and a deep null away from it. 

%convert to cartesian coordinates 
%r = 1; 


xdir = 

SdB 

. * u; %x 

coord 

normal 

to 

yz 

plane 

- amplitude 

ydir = 

SdB 

.* v; %y 

coord 

normal 

to 

xz 

plane 

- phi 

zdir = 

SdB 

. * w; % z 

coord 

normal 

to 

xy 

plane 

- theta 


figure (3) 

axis ( [0 180 0 180 40 90]) 

surf(ydir,xdir,zdir,SdB,'EdgeColornone'); %nxm 3D plot at freq 
of interest 

set(gca,'ZDir','reverse') 
camup( [ 0 1 0 ] ); 
campos( [20 15 20]) 
title(figtitle) 
xlabel('y axis') 
ylabel('x axis') %phi - n 
zlabel('z axis') %theta - m 
shading interp 

%} 

time taken = [time taken;toc]; 
end 
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beamform 12ch.m 

function beamform_12ch (Fs,binN,N,bindata) 

%use tic & toe to time beamformer cycle 
tic 

global Fpoints targetPosition targetphi targettheta time_taken; 


%DAQ 

channel assignments 




prsl 

= 

1; 

%sensor 324 pressure 



blul 

= 

2; 

%sensor 324 blu 

velocity 

Z 

axis 

grnl 

= 

3; 

%sensor 324 grn 

velocity 

X 

axis 

redl 

= 

4; 

%sensor 324 red 

velocity 

Y 

axis 

prs2 

= 

5; 

%ACO calibrated 

pressure 

mic 

blu2 

= 

6 ; 

%ACO calibrated 

pressure 

no 

velocity 

grn2 

= 

7; 

%ACO calibrated 

pressure 

no 

velocity 

red2 

= 

8; 

%ACO calibrated 

pressure 

no 

velocity 

prs3 

= 

9; 

%sensor 323 pressure 



blu3 

= 

10; 

%sensor 323 blu 

i velocity 

Z 

axis 

grn3 

= 

11; 

%sensor 323 grn 

l velocity 

X 

axis 

red3 

= 

12; 

%sensor 323 red velocity 

Y 

axis 


dummy 

dummy 

dummy 


thetainc = 1; %angle increment degrees 
phiinc = 1; %angle increment degrees 
numsens = 3; %total number of sensors 
numel = 4; %number of elements per sensor 


%calculated parameters and constants 
k = linspace(0,N-l,N); %baseline FFT bin numbers 
d = 0.172; %array element spacing meters - measured 172mm ± 
c = 340; %speed of sound meters/sec @ 20°C sea level 
maxchan = numsens * numel; %number of element channels, 6-8 
dummys here 

x = 1; %X axis designator 
y = 2; %Y axis designator 
z = 3; %Z axis designator 

%Vector Sensor Angles relative to Array Apparatus - radians 
%323 sensor 

thetabl = 0; %323 blul 0° 
phibl = pi/2; %323 blul 90° 
thetagl = pi/2; %323 grnl 90° 
phigl = 0; %323 grnl 0° 
thetarl = pi/2; %323 redl 90° 
phirl = pi/2; %323 redl 90° 

%324 sensor 

thetab3 = 0; %324 blu3 0° 
phib3 = pi/2; %324 blu3 90° 
thetag3 = pi/2; %324 grn3 90° 
phig3 = 0; %324 grn3 0° 
thetar3 = pi/2; %324 red3 90° 
phir3 = pi/2; %324 red3 90° 

%Sensor Element Unit Vectors 


2 mm 

are 
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ue(prsl, 
[omni] 

ue(blul, 
cos(thetabl) 
ue(grnl, 


) = [000]; %323 prsl dummy [000] - change in udot to 

) = [sin(thetabl)*cos(phibl) sin(thetabl)*sin(phibl) 

; %323 blul [0 0 1] 

) = [sin(thetagl)*cos(phigl) sin(thetagl)*sin (phigl) 


cos (thetagl)( 

]; %323 

grnl 

[1 0 0] 





ue(redl, 

:) = [sin(thetarl)*cos(phirl) sin(thetarl)*sin 

(phirl) 

cos(thetarl) ( 

I; %323 

redl 

[0 1 0] 





ue (prs2, 

:) = [0 

0 0]; 

%aco prs2 dummy [0 

0 

0] 

- change 

in udot to 

[omni] 








ue(blu2, 

:) = [0 

0 0]; 

%aco no vel dummy 





ue(grn2, 

:) = [0 

0 0]; 

%aco no vel dummy 





ue(red2, 

:) = [0 

0 0]; 

%aco no vel dummy 





ue(prs3, 

:) = [0 

0 0]; 

%324 prs3 dummy [0 

0 

0] 

- change 

in udot to 


[omni] 

ue(blu3, 
cos(thetab3) 
ue(grn3, 
cos(thetag3) 
ue(red3, 
cos(thetar3) 


= [sin(thetab3)*cos(phib3) sin(thetab3)*sin(phib3) 

%324 blu3 [0 0 1] 

= [sin(thetag3)*cos(phig3) sin(thetag3)*sin(phig3) 

%324 grn3 [1 0 0] 

= [sin(thetar3)*cos(phir3) sin(thetar3)*sin(phir3) 

%324 red3 [0 1 0] 


%Steered Unit Vectors 

thetas = 0:deg2rad(thetainc):pi; %polar steer angle array- radians 
phis = -pi:deg2rad(phiinc):pi; %azimuthal steer angle array upper 
hemisphere - radians 

%phis = -pi/2:deg2rad(phiinc):pi/2; %azimuthal steer angle array 
upper hemisphere - radians 

m = length(thetas) ; %total number of angles in theta direction 
n = length(phis); %total number of angles in psi direction 
p = length(k); %vector length in k direction 

u = sin(thetas') * cos(phis); %steered direction cosine x - m x n 
matrix 

v = sin(thetas') * sin(phis); %steereddirection cosine y - m x n 
matrix 

w = cos(thetas') * ones(l,n); %steereddirection cosine z - m x n 
matrix 


%element unit vector & steer angle dot product array 
%calculate one udot per element, same for all sensors, 

%use first four channels of ue matrix 

%vectorize dot products for faster multiplication in beamformer 
udot = ue(prsl:redl,x)*u(:)' + ue(prsl:redl,y)*v(:)' + 
ue(prsl:redl,z)*w(:)'; %l-4 x m*n 

%omni-directional pressure sensors - vector is "1" in all 
orientations and aspects 

udot(l,:) = ones(l,m*n); %1 x m*n 


Xc = zeros(maxchan,N); 

%hanning window to reduce sidelobes of signal 
H = hann(binN) * ones(1,maxchan); 
xwin = bindata . * H; 

%clear bindata H 
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%convert binary signal data to frequency domain 
Xw = fft(xwin); %by columns binN x maxchan 

%truncate length of signal to match transfer function length N x 
maxchan 

X = Xw(1:N,:); 

X = X.'; %change to row data maxchan x N, non conjugate transpose 
clear xwin Xw %save memory 


load '.\txfunction 323.mat' 

Hplb_BLUE = Hp1B(1:N,:); 

Hplb_GREEN = Hp1G(1:N,:); 

Hplb_RED = HplR(1:N, :) ; 

%Apply transfer functions 

Xc (1, :) = X(l, :); 

Xc(2,:) = X(2,:) . * Hplb BLUE'; %Hbl=Hpl*Hbl - pre multiply vel and 

prs Xfer funct 

Xc(3,:) = X(3,:) .* Hplb GREEN'; %Hgl=Hpl*Hgl - pre multiply vel 

and prs Xfer funct 

Xc(4,:) = X(4,:) .* Hplb RED'; %Hrl=Hpl*Hrl - pre multiply vel and 

prs Xfer funct 

Xc(9,:) = X(9,:); 

Xc(10,:) = X(10,:) . * Hplb BLUE'; %Hbl=Hpl*Hbl - pre multiply vel 

and prs Xfer funct 

Xc(ll,:) = X(ll,:) .* Hplb_GREEN'; %Hgl=Hpl*Hgl - pre multiply vel 

and prs Xfer funct 

Xc(12,:) = X(12,:) . * Hplb RED'; %Hrl=Hpl*Hrl - pre multiply vel 

and prs Xfer funct 


%apply band-pass filter 
length(Xc) 
startband = 300; 
stopband = 600; 

startbin = floor(startband*2*N/Fs); 
stopbin = floor(stopband*2*N/Fs); 

Xc2 = zeros(maxchan,N); 

Xc2 (:,startbin:stopbin) = Xc(:,startbin:stopbin); 

%Use peakfinder.m to determine source frequency/s 
[peakLocl, peakMag]=peakfinder(Xc2(1, :)) ; 

peakLoc = peakLocl.*Fs./(N*2); %convert from bin numbers to actual 
frequencies 

%choose maximum peak within earlier specified passband 
[C,I] = max(peakMag); 

Fo = peakLoc(I); 

Fpoints = [Fpoints; Fo] ; %store into global variable for appending 
onto LOFARgram 
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%+- bins around design freq to truncate - processing bandwidth, 
ignore remainder 

Fplot = floor(Fo*(N*2)/Fs) + 1; %k bin for freq Fo 

HBW = Fo/100; %+- freq to offset signal from Fo - broad for live 
signals 

halfbandwidth = floor(HBW*(N*2)/Fs) + 1; %+- number of bins to 
offset signal 

kshort = k(: , Fplot-halfbandwidth:Fplot+halfbandwidth); %get 
relevant k values 

k = kshort; %transfer relevant k values to new k vector 

clear kshort 

p = length(k); %new vector length in k direction 


%select appropriate bins for beamforming 

%Option 1: Use detected peak frequency Fo +- buffer 

Xcshort = Xc(:,Fplot-halfbandwidth:Fplot+halfbandwidth); 

%Option 2: Use continuous range of frequencies for noise reduction 
average 

%Xcshort = Xc(:,startbin:stopbin); 

Xc = Xcshort; 

clear Xcshort %save memory 
Xc = conj (Xc); 


%%initialize beamformer%% 

%%S is sum of all elements and sensors - the beamformer%% 

S = zeros(m*n,size(Xc,2)); %final beam former output - m*n x p 
array 


%Beamformer output 
for q = 9:12; 

%pull element matrix from udot array 
%Theta m * Phi n vector 

udotc = squeeze(udot((mod(q+3,4)+1); %m*n x 1 
sumcalc = udotc(:) * Xc(q,:); %m*n x p 
S = S + sumcalc; 

end 

S = sum(abs(S),2); %sum along k "amplitude", collapse to single 
value in k direction 

S = reshape(S, [m n]); %mxn thetaxphi 2D matrix 

%remove SdB points outside 90% max value for clearer phi-theta 
plots 

indices = S > max(max(S)).*0.9; 

S = indices.*S; 

%convert to dB 

SdB = 20*logl0(abs (S)); 

%Find target phi and theta by searching for maximum values 
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[SdBmaxl,11] = max(max(SdB, [ ] , 1) ) ; 

[SdBmax2,I2] = max(max(SdB, [ ],2)); 
phisdeg = rad2deg(phis); 
thetasdeg = rad2deg(thetas); 
targetphi = phisdeg(II); 
targettheta = thetasdeg(12) ; 
targetphirad = phis(II); 
targetthetarad = thetas(12); 

%Append onto global vector for final tracking plot 
targetPosition = [targetPosition; [targetphi targettheta]]; 

%Convert into coordinate system of 'Cone' subfunction 
x = sin(targetthetarad) * sin(targetphirad) * 20; 
y = cos (targetthetarad) * 20; 

z = sin (targetthetarad) * cos (targetphirad)* 20; 

Cone([0 0 0],[x y z],[0 l],10,'r',0,0,Fo,targettheta,targetphi); 


%Overhead Theta by Phi amplitude plot 
figure (5) 

%clims = [min(min(SdB)) max(max(SdB))]; 

%imagesc (rad2deg (phis) , rad2deg (thetas) , SdB, dims) ; 
imagesc(rad2deg(phis),rad2deg(thetas) , SdB) ; 

location = strcat({ ' \phi=' }, num2str(targetphi), {' \theta='}, 

num2str(targettheta) ) ; 

figtitle = strcat ({' Detected source of '}, num2str(Fo,' %4.Of' ), {' 

Hz'], location); 

title(figtitle) 
axis([-180 180 0 180]) 
xlabel( '\phi (deg)' ) 
ylabel( '\theta (deg)' ) 
shading interp 

%{ 

%3D rendering of sensor beam pattern. This allows visual 
verification 

of vector sensor's performance, which is has a maximum in the 
direction 

of source, and a deep null away from it. 

%convert to cartesian coordinates 
%r = 1; 


xdir = 

SdB 

. * u; %x coord 

normal 

to 

yz 

plane 

- amplitude 

ydir = 

SdB 

.* v; %y coord 

normal 

to 

xz 

plane 

- phi 

zdir = 

SdB 

. * w; %z coord 

normal 

to 

xy 

plane 

- theta 


figure (3) 

axis ( [ 0 180 0 180 40 90]) 

surf(ydir,xdir,zdir,SdB,'EdgeColor','none'); %nxm 3D plot at freq 
of interest 

set(gca,'ZDir','reverse') 
camup([0 1 0]) ; 
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campos([20 15 2 
title(figtitle) 
xlabel('y axis' 
ylabel('x axis' 
zlabel ( 'z axis' 
shading interp 

%} 

time taken = [time 
end 


0 ] ) 

) 

) %phi - n 
) %theta - 


taken;toe] 


SUPPORTING FUNCTIONS 


peakfinder.m 


function varargout = peakfinder(xO, sel, thresh, extrema) 

% http://www.mathworks.com/matlabcentral/fileexchange/25500 

O, 

o 

% PEAKFINDER Noise tolerant fast peak finding algorithm 

O, 

o 

% INPUTS: 

% xO - A real vector from the maxima will be found (required) 

% sel - The amount above surrounding data for a peak to be 

% identified (default = (max(xO)-min(xO))/4). Larger values 

mean 

% the algorithm is more selective in finding peaks. 

% thresh - A threshold value which peaks must be larger than to 


be 


that 


the 


maxima or smaller than to be minima, 
extrema - 1 if maxima are desired, -1 if minima are desired 
(default = maxima, 1) 

OUTPUTS: 

peakLoc - The indicies of the identified peaks in xO 
peakMag - The magnitude of the identified peaks 

[peakLoc] = peakfinder(xO) returns the indicies of local maxima 

are at least 1/4 the range of the data above surrounding data. 

[peakLoc] = peakfinder(xO,sel) returns the indicies of local maxima 
that are at least sel above surrounding data. 

[peakLoc] = peakfinder(xO,sel,thresh) returns the indicies of local 
maxima that are at least sel above surrounding data and larger 
(smaller) than thresh if you are finding maxima (minima). 

[peakLoc] = peakfinder(xO,sel,thresh,extrema) returns the maxima of 
data if extrema > 0 and the minima of the data if extrema < 0 


% [peakLoc, peakMag] = peakfinder(xO,...) returns the indicies of the 

% local maxima as well as the magnitudes of those maxima 

o, 

o 

% If called with no output the identified maxima will be plotted 
along 

% with the input data. 


Note: If repeated values are found the first is identified as the 


peak 


% Ex: 

% t = 0: .0001:10; 

% x = 12*sin(10*2*pi*t)-3*sin(.l*2*pi*t)+randn(1,numel(t)); 
% x(1250:1255) = max(x); 
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peakfinder(x) 


o, 
o 

o, 
o 

% Copyright Nathanael C. Yoder 2011 (nyoder@gmail.com) 

% Perform error checking and set defaults if not passed in 
error(nargchk(1,4,nargin, 'struct' )); 
error(nargoutchk(0,2,nargout, 'struct' )); 

s = size(xO); 

flipData = s(l) < s(2); 

lenO = numel(xO); 

if lenO ~= s(1) && lenO ~= s(2) 

error (' PEAKFINDER:Input' , 'The input data must be a vector') 
elseif isempty(xO) 

varargout = {[],[]}; 
return; 

end 

if -isreal(xO) 

warning (' PEAKFINDER:NotReal' , 'Absolute value of data will be used') 
xO = abs(xO); 

end 

if nargin <2 || isempty(sel) 
sel = (max(xO)-min(xO))/4; 
elseif -isnumeric(sel) || -isreal(sel) 

sel = (max(xO)-min(xO))/4; 
warning( 'PEAKFINDER:InvalidSel' , . . . 

'The selectivity must be a real scalar. A selectivity of %.4g 
will be used', sel) 
elseif numel(sel) > 1 

warning( 'PEAKFINDER:InvalidSel' , ... 

'The selectivity must be a scalar. The first selectivity value 
in the vector will be used.') 
sel = sel(1); 

end 

if nargin < 3 || isempty(thresh) 

thresh = [] ; 

elseif -isnumeric(thresh) || -isreal(thresh) 
thresh = []; 

warning(' PEAKFINDER:InvalidThreshold' , . . . 

'The threshold must be a real scalar. No threshold will be 

used.' ) 

elseif numel(thresh) > 1 
thresh = thresh (1); 

warning( 'PEAKFINDER:InvalidThreshold' , . . . 

'The threshold must be a scalar. The first threshold value in 
the vector will be used.') 
end 

if nargin < 4 || isempty(extrema) 
extrema = 1; 

else 

extrema = sign(extrema(1)); % Should only be 1 or -1 but make sure 
if extrema == 0 


68 


error (' PEAKFINDER:ZeroMaximaEither 1 (for maxima) or -1 (for 
minima) must be input for extrema'); 
end 

end 

xO = extrema*xO(:); % Make it so we are finding maxima regardless 
thresh = thresh*extrema; % Adjust threshold according to extrema. 
dxO = diff(xO); % Find derivative 

dxO(dxO == 0) = -eps; % This is so we find the first of repeated values 
ind = find(dxO(1:end-l) .*dx0(2 rend) < 0)+l; % Find where the derivative 
changes sign 

% Include endpoints in potential peaks and valleys 
x = [xO(1);x0(ind);x0(end)]; 
ind = [ 1;ind;lenO]; 

% x only has the peaks, valleys, and endpoints 
len = numel(x); 
minMag = min(x); 

if len > 2 % Function with peaks and valleys 

% Set initial parameters for loop 
tempMag = minMag; 
foundPeak = false; 
leftMin = minMag; 

% Deal with first point a little differently since tacked it on 
% Calculate the sign of the derivative since we taked the first 
point 

% on it does not neccessarily alternate like the rest. 
signDx = sign(diff(x(1:3))) ; 

if signDx (1) <= 0 % The first point is larger or equal to the 
second 

i i = 0; 

if signDx(1) == signDx(2) % Want alternating signs 

x(2) = []; 
ind(2) = []; 
len = len-1; 

end 

else % First point is smaller than the second 

ii = 1; 

if signDx(1) == signDx(2) % Want alternating signs 

x(1) = []; 
ind(1) = []; 
len = len-1; 

end 

end 

% Preallocate max number of maxima 
maxPeaks = ceil(len/2); 
peakLoc = zeros(maxPeaks,1); 
peakMag = zeros(maxPeaks,1); 
clnd = 1; 
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% Loop through extrema which should be peaks and then valleys 
while ii < len 

ii = ii+1; % This is a peak 

% Reset peak finding if we had a peak and the next peak is 

bigger 

% than the last or the left min was small enough to reset, 
if foundPeak 

tempMag = minMag; 
foundPeak = false; 

end 

% Make sure we don't iterate past the length of our vector 
if ii == len 

break; % We assign the last point differently out of the 

loop 

end 

% Found new peak that was lager than temp mag and selectivity 

larger 

% than the minimum to its left, 
if x(ii) > tempMag && x(ii) > leftMin + sel 
tempLoc = ii; 
tempMag = x(ii); 

end 

ii = ii+1; % Move onto the valley 
% Come down at least sel from peak 
if ~foundPeak && tempMag > sel + x(ii) 

foundPeak = true; % We have found a peak 
leftMin = x(ii); 

peakLoc(clnd) = tempLoc; % Add peak to index 
peakMag(clnd) = tempMag; 
clnd = clnd+l; 

elseif x(ii) < leftMin % New left minima 
leftMin = x(ii); 

end 

end 

% Check end point 

if x(end) > tempMag && x(end) > leftMin + sel 
peakLoc(clnd) = len; 
peakMag(clnd) = x(end); 
clnd = clnd + 1; 

elseif ~foundPeak && tempMag > minMag % Check if we still need to 
add the last point 

peakLoc(clnd) = tempLoc; 
peakMag(clnd) = tempMag; 
clnd = clnd + 1; 

end 

% Create output 

peaklnds = ind(peakLoc(1:clnd-l)) ; 
peakMags = peakMag(1:clnd-l) ; 

else % This is a monotone function where an endpoint is the only peak 
[peakMags,xlnd] = max(x); 
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if peakMags > minMag + sel 
peaklnds = ind(xlnd); 

else 

peakMags = []; 
peaklnds = []; 

end 

end 

% Apply threshold value. Since always finding maxima it will always be 
% larger than the thresh, 
if ~isempty(thresh) 

m = peakMags>thresh; 
peaklnds = peaklnds(m); 
peakMags = peakMags(m); 

end 


% Rotate data if needed 
if flipData 

peakMags = peakMags.'; 
peaklnds = peaklnds.'; 

end 


% Change sign of data if was finding minima 
if extrema < 0 

peakMags = -peakMags; 
xO = -xO; 

end 

% Plot if no output desired 
if nargout == 0 

if isempty(peaklnds) 

disp('No significant peaks found') 

else 

figure; 

plot(1:lenO,xO, ' . - ' ,peaklnds,peakMags, 'ro' , 'linewidth , 2 ); 

end 

else 

varargout = {peaklnds,peakMags}; 

end 
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cone.m 


function [Cone,EndPlatel,EndPlate2] = 

Cone(XI,X2,R,n,cylcolor,closed,lines,Fo,targettheta,targetphi) 

O, 

O 

% This function constructs a cylinder connecting two center points 

Q, 

O 

% Usage : 

% [Cone,EndPlatel, EndPlate2] = Cone(XI,X2,R,n,cyl_color,closed,lines) 

o, 

o 

% Cone-Handle of the cone 

% EndPlatel-Handle of the Starting End plate 

% EndPlate2-Handle of the Ending End plate 

% XI and X2 are the 3x1 vectors of the two points 

% R is the radius of the cylinder/cone R(l) = start radius, R(2) = end 
radius 

% n is the no. of elements on the cylinder circumference (more--> 
refined) 

% cyl color is the color definition like 'r','b',[0.52 0.52 0.52] 

% closed=l for closed cylinder or 0 for hollow open cylinder 
% lines=l for displaying the line segments on the cylinder 0 for only 
% surface 

o, 

o 

% Typical Inputs 
% Xl=[10 10 10]; 

% X2=[35 20 40]; 

% r= [ 1 5]; 

% n=20; 

% cyl_color='b'; 

% closed=l; 

o, 

o 

% NOTE: There is a MATLAB function "cylinder" to revolve a curve about 
an 

% axis. This "Cylinder" provides more customization like direction and 
etc 


% Calculating the length of the Cone 
length_cyl=norm(X2-X1) ; 

% Creating 2 circles in the YZ plane 

t=linspace(0,2*pi,n) 

xa2=R(1)*cos(t); 

xa3=R(1)*sin(t); 

xb2=R(2)*cos (t); 

xb3=R(2)*sin (t); 

% Creating the points in the X-Direction 
xl=[0 length_cyl]; 

% Creating (Extruding) the cylinder points in the X-Directions 
xxl=repmat(xl,length(xa2) , 1) ; 
xx2=[xa2 xb2]; %xx2=repmat(x2,1,2); 
xx3=[xa3 xb3]; %xx3=repmat(x3,1,2) ; 
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% Drawing two filled cirlces to close the cylinder 
if closed==l 
hold on 

EndPlatel=fill3(xxl(:,1),xx2(:,1),xx3(:,1), ' r ' ) ; 
EndPlate2=fill3(xxl(:,2),xx2(:,2) , xx3(: , 2) , ' r ' ) ; 

end 


% Plotting the cylinder along 
starting 
% from Origin 
figure (2) 

Cone=mesh(xxl,xx2,xx3); 

D = 20; 

axis([-D D -D D -D D]) 
%xlabel('y axis') 

%ylabel('x axis') 

%zlabel('z axis') 

camup([0 0 1] ) 

campos([70 15 20]) 

location = strcat ({ ' \phi=' }, 

num2str(targettheta)); 

figtitle = strcat ({' Detected 

Hz'}, location) ; 

title(figtitle) 


the X-Direction with required length 


num2str(targetphi), {' \theta='}, 

source of num2str(Fo,' %4.Of' ), {' 


% Defining Unit vector along the X-direction 
unit_Vx=[l 0 0]; 


% Calulating the angle between the x direction and the required 
direction 

% of Cone through dot product 

angle XlX2=acos( dot( unit_Vx,(X2-X1) )/( normfunit Vx)*norm(X2-X1)) 

)*1807pi; 


% Finding the axis of rotation (single rotation) to roate the Cone in 
% X-direction to the required arbitrary direction through cross product 
axis_rot=cross([1 0 0],(X2-X1) ); 

% Rotating the plotted Cone and the end plate circles to the required 
% angles 

if angle X1X2~=0 % Rotation is not needed if required direction is 
along X 

rotate(Cone,axis_rot,angle_XlX2,[0 0 0]) 
if closed==l 

rotate(EndPlatel, axis rot,angle X1X2, [0 0 0]) 
rotate(EndPlate2,axis_rot,angle X1X2,[0 0 0]) 

end 

end 


% Till now Cone has only been aligned with the required direction, but 
% position starts from the origin, so it will now be shifted to the 
right 

% position 
if closed==l 

set(EndPlatel, 'XData ,get(EndPlatel, XData' )tXl(1)) 
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set(EndPlatel, 'YData ,get(EndPlatel,' YData' )+X1(2)) 
set(EndPlatel, 'ZData' ,get(EndPlatel,' ZData' )+X1(3)) 

set(EndPlate2, 'XData ,get(EndPlate2, XData' )+X1(1)) 
set(EndPlate2, 'YData' ,get(EndPlate2,' YData' )+X1(2)) 
set(EndPlate2 ZData ,get(EndPlate2, ZData' )+X1(3)) 

end 

set(Cone, 'XData 1 ,get(Cone,' XData' )+X1(1)) 
set(Cone, 'YData' , get(Cone, 'YData' )+X1(2)) 
set(Cone, 'ZData ',get(Cone, 'ZData' )+X1(3)) 

% Setting the color to the Cone and the end plates 

set(Cone, 'AmbientStrength' ,1, 'FaceColor' ,cyl color, 'FaceLighting' , 'gour 
aud' ); %,'EdgeColor','none') 
if closed==l 

set([EndPlatel 

EndPlate2], 'AmbientStrength ',1,'FaceColor', cyl color, 'FaceLighting' , 'go 

uraud EdgeColor','none') 

else 

EndPlatel= [ ]; 

EndPlate2=[]; 

end 

% If lines are not needed making it disapear 
if lines==0 

set(Cone, 'EdgeAlpha' ,0) 

end 

%shading faceted % faceted flat interp; 

%camlight; 

%light; 

%lighting gouraud; %flat gouraud phong none 
material dull; %shiny dull metal 
%surcolormap(bone) 

%camlight headlight; 

%light('Style', 'local', 'Position', [720 0 500]); 
light ('Style', 'local','Position', [0 480 500]); 
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o\° o\° 


D. 


MISCELLANEOUS 


CPA rootfinder.m 

%Two methods for solving the horizontal distance at CPA. 


n = 

1; 

%nth frequency maxima 

c = 

343; 

%speed of sound 

h = 

200; 

%target altitude 

d = 

0.95; 

%sensor height 

fn ^ 

= 182; 

%frequency 


Exact solution by root-finding. Take the positive root for obvious 
reasons. 

f = @(x)sqrt((h+d). A 2 + x. A 2)-sqrt((h-d). A 2 + x. A 2) - n*c/fn; 

r = fzero(f,10) 

%Approximation when r << h+d 
r2 = h*sqrt(2-(n*c)/(fn*d)) 
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